Time Dependent Behavior of 
Granular Material in a Vibrating Box 



Jysoo Lee 

Benjamin Levich Institute and Department of Physics 
City College of the City University of New York 
New York, NY 10031 

February 1, 2008 

Abstract 

Using numerical and analytic methods, we study the time depen- 
dent behavior of granular material in a vibrating box. We find, by 
molecular dynamics simulation, that the temporal fluctuations of the 
pressure and the height expansion scale in Af, where A (/) is the am- 
plitude (frequency) of the vibration. On the other hand, the fluctua- 
tions of the velocity and the granular temperature do not scale in any 
simple combination of A and /. Using the kinetic theory of Haff, we 
study the temporal behaviors of the hydrodynamic quantities by per- 
turbing about their time averaged values in the quasi-incompressible 
limit. The results of the kinetic theory disagree with the numeri- 
cal simulations. The kinetic theory predicts that the whole material 
oscillates roughly as a single block. However, the numerical simula- 
tions show that the region of active particle movement is localized and 
moves with time, behavior very similar to the propagation of a sound 
wave. 
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1 Introduction 



Systems of granular particles (e.g. sand) exhibit many interesting phenom- 
ena, such as segregation under vibration or shear, density waves in the out- 
flow through a hopper and a tube, and the formation of heaps and convection 
cells under vibration [1-5]. These phenomena are consequences of the unusual 
dynamical response of the systems, and are for the most part still poorly un- 
derstood. 

We focus on the vertical vibration of a box containing granular parti- 
cles. There are many interesting phenomena associated with this system, 
such as convection cells [6-11], heap formation [12-17], sub-harmonic insta- 
bility f[8 |, surface waves |19|, ^(J and even turbulent flows [ZT|. The basis 
for understanding these diverse phenomena is the state of granular media 
under vibration. The state is characterized by the hydrodynamic fields of 
the system, such as the density, velocity and granular temperature fields. 

There have been several studies on the state of granular particles under 
vibration. Thomas et al. studied the system in three dimensions, mainly 
focusing on the behavior of shallow beds f22|. Clement and Rajchenbach 
experimentally measured the density, velocity and temperature fields of a 
two-dimensional vertical packing of beads ||2~3| . They found that the tem- 
perature increases monotonically with the distance from the bottom plate. 
The same system was studied by molecular dynamics (MD) simulation with 
similar results ||24j| . In a series of simulations and experiments, Luding et 



2 



al. studied the behavior of the one and two-dimensional systems [25-27]. 
They found that the height expansion, which is the rise of the center of mass 
due to the vibration, scales in the variable x = Af. Here, A and / are 
the amplitude and the frequency of the vibration. Warr et al. experimen- 
tally confirmed the scaling, and they also gave an argument for its origin 
. In recent MD simulations of the three-dimensional system, Lan and 



Rosato measured the density and temperature fields p9| . They compared 
the results with the theoretical predictions by Richman and Martin |3(| , and 
found good agreement. Also, an approximate theory was developed for the 
system in one dimension, which agrees with simulations in the weak and the 



strong dissipative regimes 31 



In a previous paper, we studied the time averaged behavior of the two- 
dimensional system using numerical and analytic methods f3"2|| . Using MD 
simulation, we found that the time averaged value of not only the expansion 
but also the density and the granular temperature fields scale in x. We 
also used the kinetic theory of Haff |33| to determine the time averaged 
hydrodynamic fields in the quasi- incompressible limit. The results are, in 
general, consistent with the numerical data, and in particular show scaling 
behavior in the variable x. We found that the origin of the scaling can be 
understood within the framework of the theory. 

In the present paper, we extend our study to the time dependent be- 
havior of the system, whose understanding is not only essential in studying 
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various time dependent phenomena, but also necessary to understand the 
mechanisms of certain steady state phenomena. For example, many of the 
arguments for the mechanism of the convection involve the variations of cer- 
tain hydrodynamic quantities over a vibration cycle [0, §, |H], TT|, [16| . 



We find, by MD simulation, that the temporal fluctuations of the pressure 
and the height expansion scale in x, while the velocity and the granular tem- 
perature do not scale in any simple combination of A and /. We also study 
the system using the kinetic theory of Haff, where the temporal behavior 
is studied by perturbing about the time average in the quasi-incompressible 
limit. The results of the kinetic theory disagree with the numerical data. 
The kinetic theory predicts that the whole system of particles moves as one 
effective "block." The numerical simulations, on the other hand, show that 
the region of active particle movement is localized, and moves with time. 
The presence of the wave is partly responsible for the discrepancy. 

The paper is organized as follows. In Sec. [2], we specify the interaction 
of the particles used in the MD simulations. We then present the temporal 
fluctuations of the expansion and the fields obtained by the simulations. 
Analytic results will be discussed in Sec. ||. The continuum equations for 
granular material will be given, and perturbation equations are derived. We 
present the solution of the equations, and compare with the numerical data. 
In Sec. (|, we check the assumptions used in obtaining the solution. We also 
discuss various properties of the waves. Conclusions are given in Sec. |5|. 
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2 Numerical Simulation 



We start by describing the interactions used in the MD simulations. The 
simulations are done in two dimensions with disk shaped particles. The 
interaction between the particles is that of Cundall and Strack |34}1 , which 



allows the particles to rotate as well as translate. Particles interact only by 
contact, and the force between two such particles i and j is the following. 
Let the coordinate of the center of particle % (j) be Ri (Rj), and r '— Ri — Rj. 
We use a new coordinate system defined by two vectors n (normal) and s 
(shear). Here, n = f/\f\, and s is obtained by rotating n clockwise by 7r/2. 
The normal component F?_+i of the force acting on particle i from particle j 
is 

■tfjlrf = K{&i + aj - \f\) - ^ n m e (v ■ h), (1) 

where a« (aj) is the radius of particle % (j), and v = df/dt. The first term 
is the linear elastic force, where k n is the elastic constant of the material. 
The constant 7„ of the second term is the friction coefficient of the velocity 
dependent damping force, and m e is the effective mass, mim,j/{mi + rrij). 
The shear component Fj_^ is given by 

Fj-M = -sign(5 S ) min(^|^|,/i|^|). (2) 

The term represents static friction, which requires a finite amount of force 
({iFj 1 ^) to break a contact. Here, \x is the friction coefficient, 5s the total 
shear displacement during a contact, and k s the elastic constant of a virtual 



tangential spring. 

The shear force also affects the rotation of the particles. The torque 
acting on particle i due to particle j is 

Tj-*i = f c xs FJ^, (3) 

where r c is the vector from the center of particle i to the point where particles 
i and j overlap. Since the particles used in the simulations are very stiff (large 
k n ), the area of the overlap is very small. It is thus a good approximation to 
use — a;ji as f c . 

A particle can also interact with a wall. The force and torque on particle 
i, in contact with a wall, are given by (jl]) - © with aj = and m e = m^. A 
wall is assumed to be rigid, i.e., it is not moved by collisions with particles. 
Also, the system is under a gravitational field g. A more detailed explanation 
of the interaction is given elsewhere 



The movements of the particles are calculated using a fifth order predictor- 
corrector method. We use two Verlet tables. One is a usual table with fi- 
nite skin thickness. The other table is a list of pairs of actually interacting 
particles, which is needed to calculate the shear force. The interaction pa- 
rameters used in this study are fixed as follows, unless otherwise specified: 
k n = 10 5 ,k s = 10 5 ,7 n = 2 x 10 2 and /i = 0.2. The timestep is taken to 
be 5 x 10 -6 . This small timestep is necessary for the large elastic constant 
used in the simulations. For too small values of the elastic constant, the 
system loses the character of a system of distinct particles, and behaves like 
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a viscous material. In order to avoid artifacts of a monodisperse system 
(e.g., hexagonal packing), we choose the radius of the particles from a Gaus- 
sian distribution with the mean 0.1 and the width 0.02. The density of the 
particles is 0.1. Throughout this paper, CGS units are implied. 

We put the particles in a two-dimensional rectangular box. The box 
consists of two horizontal (top and bottom) plates which oscillate sinusoidally 
along the vertical direction with given amplitude A and frequency /. The 
separation between the two plates H is chosen to be much larger (10 5 times) 
than the average radius of the particles, so the particles do not interact 
with the top plate for all cases studied here. We apply a periodic boundary 
condition in the horizontal direction. The width of the box is W — 1. We 
also try different values of W, and find no essential difference in the following 
results. 

We start the simulation by inserting the particles at random positions 
in the box. We let them fall by gravity and wait while they lose energy by 
collisions. We wait for 10 5 iterations for the particles to relax, and during 
this period we keep the plates fixed. The typical velocity at the end of the 
relaxation is of order 10" 2 . After the relaxation, we vibrate the plates for 50 
cycles before taking measurements in order to eliminate any transient effect. 
Measurements are made during the next 200 cycles. 

We measure hydrodynamic quantities — density, velocity and granular 
temperature — which characterize the state of the system. The most detailed 
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information is contained in the time series of their fields (e.g., density field), 
which will be discussed later. Here, we want to start with something simple 
and representative of the system. 

The center of mass of the particles can be loosely related to the density. 
Let y(t) be the vertical coordinate of the center of mass at time t. The mean 
density is related to the mean interparticle distance, which is also related to 
y(t). Since y(t) is a scalar which also can be easily measured, we study it 
as a representative of the density. In the same spirit, we study the space 
averaged vertical velocity V y (t) and granular temperature r(i) instead of the 
complete local fields. We define the spatial average of A(x, y, t) as 

f f A(x,y,t)p(x,y,t)dxdy 
{A{t)) = oxytdxdy ~ (4) 

where p(x,y,t) is the density field at position (x,y). It thus follows that 
V y (t) = (v y (t)) and r(t) = (T(t)), where v y (x,y,t) and T(x,y,t) is the 
vertical velocity and the granular temperature field, respectively. Another 
quantity of interest is the pressure field p(x,y,t). For information on this 
quantity, we will study the total pressure at the bottom 

_ fp(x, y = 0, t)p(x, y = 0, t)dx 
M) ~ Jp(x,y = 0,t)dx ■ [b) 

We studied the time averaged behavior of the system in our previous 



paper [j32| . Here, we study the temporal behavior, especially focusing on 
the variations of the fields within a vibration cycle. We first measure the 
temporal fluctuation, which is defined as the standard deviation of a temporal 



sequence. Let (A) t be the time average of A(t); then, we use 

l/cxp = (ycxp)t, 5l/ exp = yj (l/cxp)i — (?/cxp)t 



% = (Vy) t , 

r = (r) t , 





- ( V y)t 








~ {Po)l 



Po = (Po)t, Sp = (p 2 ) t - (p )l (6) 

where the expansion y e -x_ P {t) is defined as the difference between y(t) during 
and before the vibration. We measure these quantities at every 1/100 of a 
period for 200 cycles. 

In our previous paper, it was shown that y exp scales in Af, in agreement 
with earlier simulations and experiments [25-28] . We find that the rotation of 
the particles included in the present simulation does not change this scaling, 
but does significantly decrease the value of y exp The decrease is probably 
due to the fact that the average translational energy becomes smaller, since 
some of the energy is transferred to the rotation. We then consider the 
fluctuation of the expansion 5y exp . In Fig. 1(a), we show 5y exp for several 
values of A and /. The quality of the scaling is not very good, but the data 
is still consistent with Af scaling, especially without the persistent deviation 
at low A part of the / = 20 data. This deviation is also present in the scaling 

Of ?/exp- 

The behavior of r(t) is very similar. It was shown that r scales in Af 



32]. Again, we find that the rotation does not change the scaling of r, 
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but does change its value. The situation becomes a little different when we 
consider the fluctuation St as shown in Fig. 1(b). It is clear from the figure 
that St does not scale in Af. In fact, it does not scale in any of the simple 
combinations of A and / we have tried. 

In the previous paper, the behaviors of V y (t) and p (t) were not discussed 
in detail, since their time averaged values are trivial. Since the system is in 
a steady state, V y is zero, which is also confirmed by the simulations. Since 
the pressure p (t) is caused by the weight of the particles, one might guess p a 
is simply the total weight of the particles divided by the area of the bottom 
(W = 1). We find that p Q is indeed a constant independent of A and /, whose 
value is consistent with the total weight. On the other hand, the behavior 
of their fluctuations is far from trivial. The fluctuation of the velocity SV y 
for several values of A and / is shown in Fig. 1(c). It is apparent that SV y 
does not scale in Af. Also, SV y does not scale in any simple combination of 
A and /. The data for the fluctuation of the pressure Sp Q is consistent with 
scaling in Af as shown in Fig. 1(d). The quality of collapse is again not very 
good, especially for the low A part of the / = 20 data. 

The behavior of the time averaged quantities is either trivial (V y and 
po), or scales in Af (y exp and f). The reason for the trivial behavior has 
been discussed, and the scaling in Af can be understood from a kinetic 
theory of granular particles |32| . The behavior of the fluctuations is, however, 
not easy to understand. For example, one might naively expect that SV y 
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behaves as Af, the velocity fluctuation of the bottom plate; but the numerical 
data suggests this is not so. Also, one might guess 5p Q behaves as Af 2 , 
the variation of the effective gravity; but the numerical simulations suggest 
scaling in Af. Thus, the behaviors of the temporal fluctuations seem to be 
inconsistent with an intuitive picture. In the next section, we discuss an 
attempt to understand this behavior. 

3 Kinetic Theory 

In this section, we study the time dependent behavior of the system using 
a kinetic theory of granular material. We use the formalism by Haff [0^1 . 
which was successfully applied to the time averaged behavior of the system 
[ j32|] . For more details on the formalism and other kinetic theories of granular 
material, see Ref. |3!| and references therein. 

Haff 's formulation consists of equations of motion for mass, momentum 
and energy conservation. The mass conservation equation is 

+ V ■ (pv) = 0, (7) 

where p and v are the density and the velocity fields, respectively. Next is 
the z-th component of the momentum conservation equation, 
d -> d -> dv ■ Ov ■ 

+ p(v- vfo = ^l-p + A(v • *)] + + ^)] + m, (8) 

where summation over index j is implied. The coefficients A and rj are vis- 
cosities which will be determined later. Also, p is the internal pressure, and 
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gi is the i-th component of the gravitational field. Although (El) resembles 
the Navier-Stokes equation, the coefficients as well as the internal pressure 
are now functions of the fields instead of being constant. The last of the 
equations of motion is energy conservation, 



OXj OX; OX a Ox,; 



Here, T is the granular temperature field, K is the "thermal conductivity," / 
is the rate of the dissipation due to inelastic collisions, and summations over 
indices i and j are implied. Although the form of (|Sp is somewhat different 
from that of the Navier-Stokes equations, the equation can still be easily 
understood. The left hand side of @ is simply the material derivative of 
the total kinetic energy, where the total kinetic energy is divided into the 
convective part (involving v) and the fluctuating part (involving T). On the 
right hand side of the equation, the first three lines are simply the rate of work 
done by the internal pressure, viscosity and gravity, respectively. The term 
involving K is the rate of energy transported by "thermal conduction." The 
dissipation term /, which is a consequence of the inelasticity of the particles, 
is responsible for many of the unique properties of granular material. 
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We now discuss the coefficients which are yet to be determined. Deriva- 
tion of the relations of these coefficients to the fields is based on intuitive 
arguments p3 |. Also, the derivation assumes that the density is not sig- 



nificantly smaller than the close-packed density, i.e., the system is almost 
incompressible. The relation for the internal pressure is 

T 

p = tdp-, (10) 

where t is an undetermined constant, and d is the average diameter of the 
particles. The variable s, which is roughly the gap between the particles, is 
related to the density by 

where m is the average mass of the particles. Then, the viscosity 77 is given 
as 



v = qd 2 P ^f, (12) 

where q is an undetermined constant. In a similar way, the thermal conduc- 
tivity is found to be 



K = rd 2 ^. (13) 
s 

Here again, r is an undetermined constant. Finally, the rate of dissipation is 

T 3/2 

I = 19 , (14) 

s 

where 7 is an undetermined constant. The viscosity A is left undetermined, 
due to the fact that, in the range where these relations are valid, the term 
containing A is negligible and is dropped from the calculation. 
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We impose two constraints in order to make the equations analytically 
tractable. The first is the horizontal periodic boundary condition. Due to 
the boundary condition, there are no significant variations of the fields along 
the horizontal direction. Thus, we only have to deal with a one- dimensional 
equation, instead of a two or three-dimensional one. The other constraint is 
incompressibility, which is a little tricky. Incompressibility implies, strictly 
speaking, that the density p is constant. Due to the relation between p and s 
([ITD) s also has to be constant. Here, we are interested in the situation where 
s is much smaller than d, but still non-zero. In such the variation of 

the density can be ignored, but not the variation of a variable that depends 
directly on s. We call this condition quasi-incompressibility. 

Under these conditions, we solved the equations for the time averaged 
fields, 

T®(y) ~ B 2 v 2 w -^— exp(-2y/£) 

Vo-y 

s (0) (y) * ^ £ 2 vl —^- 2 exp(-2y/£) 

g [Vo - vr 

vf\y) = 

P {0 \y) = Pog(y -y). (15) 

Here, v w = 2nAf the maximum velocity of the bottom, y the height of the 



free surface, and i — v/r/7 d the dissipation length. Also, p Q is the density 
of the maximum packing, and 

r>2 __ (1 + e w ) 2 

l-el-(2rd/a£)(l-£/2y o y { ) 
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where e w is the coefficient of restitution of collisions between a particle and 
a wall . 

We study the time dependent behavior of a quantity by perturbing it from 
its time averaged value. We assume that the perturbation term oscillates with 
/ — the frequency of the vibration. In general, the assumption is not valid, 
since one has to consider all the modes with different frequencies. However, 
when the amplitude of the vibration is small enough, in many cases, the mode 
with frequency / dominates the time dependent behavior. We expect there 
is a range of A in which the assumption is valid, which will be determined 
later by the numerical simulations. We thus use 



T(y,t) 


= T^(y)+T^(y) 


■ exp(iut) 




s(y,t) 


= S M(y) +S W(y). 


exp(icut) 




v y {y,t) 


= vf\y) + 4\y)- 


exp(iuot) 




p(v,t) 


= p {0 \y)+pU(y)- 


exp(iuit). 


(17) 



Substituting (|T7| ) into the mass conservation condition (|7|) and using the 
quasi- incompressibility condition, we obtain dv^\y)/dy = 0. Here and in 
the rest of the calculation, we consider terms up to the first order in the 
expansion. Since v y (y = 0,t) should be the velocity of the bottom plate, 

V W(y) = iv w , (18) 

which is a consequence of quasi-incompressibility and the one- dimensional 
nature of the system. Also, momentum conservation equation (|8|), combined 
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with (|TT|), becomes 

p il \y)=p Auj 2 (y -y), (19) 

which can be easily understood. The pressure at the bottom is proportional 
to the total weight of the particles. Thus, the fluctuation of p (t) is the 
total mass p y times the fluctuation of the effective gravity Auj 2 . Finally, 
we consider the energy conservation equation (|9|). Combined with (0), it 
becomes 

y dy 
= r A^ ip (o)d ui i )+pi i ) d u(0) 
t dy dy dy 

_ ^ ([/ (iyo) + f/ (oyi) ); (20) 

where we introduce the variable U(y, t) = \jT(y, t). One has to solve ( p0|) 
for T(y,t). We can not, unfortunately, get an analytic solution of the result- 
ing nonlinear differential equation. Since s(y,t) has to be calculated from 
the relation fllQ|) between T(y,t), s(y,t) and p(y,t), we also can not get an 
expression for s^\y). 

We compare the results from the kinetic theory with the numerical simu- 
lations (Fig. 1). First, we consider v y (y,t). Since the temporal fluctuation of 
v y (y,t) is proportional to v^(y), the kinetic theory predicts SV y ~ Af. The 
data from the simulations, however, is inconsistent with this scaling. The 
situation is similar for p(y,t). The kinetic theory predicts 5p Q ~ Af 2 , while 
the numerical data scales in Af. We do not have analytic expressions for 
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T(y,t) and s(y,t) to compare with the simulation data. 

The failure of the kinetic theory, when applied to the time dependent 
behavior, is in sharp contrast with its success in studying the time averaged 
behavior. The kinetic theory correctly predicts the scaling behavior of all the 
time averaged hydrodynamic quantities. We suspect the failure is associated 
with the breakdown of a key assumption(s) used in the theory. Three key 
assumptions, besides the ones employed in the formulation of the kinetic 
theory, are made to obtain a simple analytic solution for the time dependent 
behavior. The first is quasi- incompressibility, which is shown to be valid for 



Af < 1.5 from the study of the time averaged behavior f3^| . However, the 
scaling behavior of the time averaged quantities remains unchanged even in 
the compressible regime. The second assumption is that the time dependence 
of a quantity is a sinusoidal oscillation with frequency /, which we expect 
to be valid for small A. The last assumption is that the time dependent 
term in flTTD is much smaller than its time averaged value in order for the 
perturbation to be valid. In the next section, we check the validity of these 
assumptions by comparing to the numerical simulations. 

4 Validity of Assumptions 

First, we check the assumption that the mode with the driving frequency 
dominates the time dependent behavior of the hydrodynamic quantities. We 
obtain a time series of y e xp(t) by measuring it at every 1/100 of a period for 
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200 cycles. We then calculate the power spectrum of y e xp(t)—y ex p using a FFT 
routine in the NAG library (c06gbf). The results with / = 100 are shown in 
Fig. 2. The mode with / = 100 is dominant for A = 0.01 (Fig. 2(a)). When 
A is increased further to 0.03, however, the / ~ 20 mode becomes dominant 
(Fig. 2(b)). The results with / = 20 are entirely similar. The / = 20 mode 
loses its dominance when A is increased to about 0.5. In both cases, the 
mode with frequency / is dominant until T ~ 10, where T = Au 2 /g. 

The measurements of V y (t) andp D (t) also support the above observations. 
In Fig. 3(a), we show the variation of V y (t) in one cycle, where / = 100 and 
the data are averaged over 200 periods. The curve with A = 0.01 is nearly 
sinusoidal, while clear deviation is seen for A = 0.03. An important point 
to note is the behavior of the maximum value of V y (t). The kinetic theory 
predicts, in (|T8|), the maximum value to be proportional to A, which is clearly 
not consistent with the data. The measured value is quite smaller than what 
is predicted. For example, for A = 0.03, the predicted value is 6ti, while the 
measured value is about 2.4. In Fig. 3(b), we show p a (t) in one period, where 
again / = 100 and the data is averaged over 200 cycles. The curve seems to 
deviate from the sinusoidal even for small A, and it is difficult to determine 
whether the mode with / = 100 dominates. Again, the predicted behavior 
of the maximum value of p (t) is not consistent with the measurement. The 
maximum value of p {t) is predicted to increase linearly with A, as in (|19|), 
which is clearly not consistent with the data. The observed maximum value 
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(~ 400) is quite smaller than the predicted value (~ 4, 000). 

The mode with the driving frequency dominates the time dependent be- 
havior for small values of A, where a rough criterion for the dominance is 
T <C 10. Also, it was shown that quasi-incompressibility is valid for Af < 1.5 
[j3^| . The validity of the linear perturbation approximation (|i~7f ) is a little 
tricky. We require that the perturbation terms are smaller than the time 
averaged terms, which is valid for the expansion and the temperature. The 
two terms are, however, comparable for the pressure even at small value of 
A = 0.05. The consequence of the large pressure fluctuation on the validity 
of the perturbation is unclear. For large A, all of the above assumptions 
are not valid, which complicates the analysis of the system. For example, in 
order to study the time dependent behavior, one has to consider additional 
modes with different frequencies. 

The surprise is that the predictions for the maximum values of the vertical 
velocity and the pressure are not correct even when all the assumptions seem 
to be valid (e.g., A = 0.01 and / = 100). We inspect again the predictions of 
the kinetic theory. As given in (0), the velocity field is uniform, and its value 
is v w the velocity of the bottom. Therefore, the solution of the perturbation 
expansion suggests that the whole system of particles is moving as a single 
"block" attached to the bottom. The spatial and temporal variations of the 
other fields do not change the single block picture, but rather describe the 
structure of the block. A consequence of the picture is that the pressure at 
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the bottom is proportional to the effective gravity which reaches its maximum 
(r + l)g at phase 37r/2 of the vibration, which is exactly (|T9|). 

However, the measurements of V y (t) and p (t) are not consistent with 
the picture. The measured maximum value of V y (t) is much smaller than 
what is predicted, which suggests that only a small fraction of the particles 
move together at a given time. Also, the fact that the measured maximum 
value of p {t) is smaller than the prediction also supports this observation. 
Furthermore, the phase at which p (t) reaches the maximum is about 0, in 
contrast to 3n/2 suggested by the single block picture. The discrepancy can 
be understood as follows. Since the maximum acceleration of the bottom 
is larger than that of gravity, particles initially lying on the bottom will be 
"launched" at a certain phase of the vibration. The maximum pressure at 
the bottom will occur when most of the launched particles come back and 
collide with the bottom, which occurs around = 0. 

The direct evidence against the single block picture is the time evolution 
of the whole velocity field shown in Fig. 4(a). Here, we use / = 100, A = 0.01, 
at which the assumptions used to derive the predictions of the kinetic theory 
seem to be valid. It is clear from the figure that the region of significant mo- 
tion is localized, and travels like a wave. The propagation of the disturbance 
seems to be very similar to that of sound waves in a gas. Also, the maximum 
velocity is about 6, close to the prediction 2tt. The localization of the particle 
motion can also be seen in the time evolution of the granular temperature 
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field shown in Fig. 4(b). Again, the region of high temperature is localized, 
and travels upwards |37|, |38|. Furthermore, the location of the high temper- 
ature region coincides with that of the large velocity region. The density 
field, on the other hand, does not vary significantly as shown in Fig. 4(c), 
which agrees with the previous experiment |23| . It is clear that the presence 
of the "waves" changes the behaviors of the fields, and possibly their scaling 
properties. The absence of the waves in the single block solution is, at least, 
partly responsible for the failure of the kinetic theory. The absence is due to 
quasi- incompressibility. In fact, it is easy to derive the single block picture 
only from the quasi-incompressibility condition and one-dimensional nature 
of the equation. It is thus necessary to consider the general case of the kinetic 
theory of a compressible gas, which unfortunately is quite complicated. 

We want to finish this section by discussing some properties of the waves. 
We first consider the motion of the maximum disturbance. In Fig. 5(a), we 
show the phase 4> mauX (y) at which v y (y,t) reaches a maximum with / = 100 
and several values of A. In other words, we plot the position of the maximum 
velocity in the y — <f> plane. The velocity of the wave, inversely proportional 
to the slope of ma x(?/) curve, is a bit small. The time needed for the wave to 
propagate from the bottom to the top of the pile is of the order of the period of 
the vibration. The simulations with / = 50 show that, for the same values of 
A, the velocity of the wave does not change significantly, suggesting that there 
is a fixed time scale for the wave propagation. Also, it can be seen from the 
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figure that the velocity decreases when either y or A increases with the other 
parameters fixed. The decrease probably results from the decrease of the 
collision frequency between particles, due to the decrease of the density. Also, 
the location of the maximum granular temperature in the y— <fi plane is shown 
in Fig. 5(b), where one can see the close correlation with Fig. 5(a). In fact, the 
maximum temperature always occurs just above the maximum velocity at a 
given phase, which probably is the point of the largest velocity gradient. We 
now discuss the values of the maximum disturbances. In Fig. 6(a), we show 
the maximum value of v y (y,t) for given y — v y (y, (fi max (y))- The maximum 
value decreases with y roughly as an exponential, and the rate of the decrease 
is larger for larger A with fixed /. The decrease is due to two causes: (1) 
some of the vertical velocity component is transferred to the horizontal one 
by interparticle collisions; (2) the kinetic energy is lost by inelastic collisions. 
The results with / = 20 are essentially the same. Also, the behavior of the 
value of the maximum temperature T(y,(f> max (y)), as shown in Fig. 6(b), is 
very similar to that of v y (y, (f) max (t)). It decreases roughly as an exponential, 
and it decays faster with larger A with fixed /, where the origin of the 
decrease is the same as the velocity field. 

5 Conclusion 

We have studied the temporal behavior of granular material in a vibrating 
box. We find that the temporal fluctuations of y cxp {t) and p (t) scale in 
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Af, while no scaling is found for the fluctuations of r(£) and V y (t). We 
study the behavior using the kinetic theory of Haff, where we perturb the 
hydrodynamic quantities from their time averaged values. The results of the 
kinetic theory are not consistent with the numerical data. We argue that 
the failure of the theory is, at least, partly due to the waves found in the 
simulations. 

We discuss some possible ways to study the time dependent behavior of 
the system. Since the main problem of the present theory, as discussed above, 
is the condition of quasi-incompressibility, it sounds reasonable to study the 
system of equations in a fully compressible regime. However, the system of 
equations becomes too complicated. The three conservation equations 
are written in general form, and do not need any modification. The relations 
of p,rj,K,I to the fields (]T0|)-(|T4|) as well as the boundary conditions have 
to be modified. It is not the modifications themselves, but the complexity 
of the resulting equations, that makes an analytic solution too difficult to 
obtain. However, we can still gain some information about the system by a 
perturbative or an approximate method, as well as the numerical solution of 
the kinetic equations. 

We thank Joel Koplik for many useful discussions and comments on the 
manuscript, and Michael Tanksley for critical reading of the manuscript. 
This work is supported in part by the Department of Energy under grant 
DE-FG02-93-ER14327. 
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Figure Captions 

Fig. 1: Scaling behaviors of the hydrodynamic quantities. Each datum is 
averaged over at least 3 samples, where 20000 measurements are made 
in a sample, (a) The fluctuation of the expansion 5y cxp seems to scale 
in Af, where (b) the fluctuation of the temperature St does not scale, 
(c) The fluctuation of the velocity 8V y does not scale, where (d) the 
fluctuation of the pressure 8p seems to scale in Af. 

Fig. 2: Power spectrum of y e x P (t) with / = 100 and (a) A = 0.01, (b) 
A = 0.03. The expansion is measured at every 10~ 4 second for 200 
cycles. The mode with / = 100 is dominant at A = 0.01, but loses its 
dominance at A = 0.03. 

Fig. 3: Time evolution of (a) V y (t), (b) p (t) in one cycle, where / = 100 
and the data are averaged over 200 cycles. Deviation from a sinusoidal 
is apparent for A = 0.03. 

Fig. 4: Time evolution of (a) v y (y,t), (b) T(y,t) and (c) p(y,t) fields, where 
/ = 100, A = 0.01 and the data are averaged over 200 cycles. The 
density field is normalized to be the volume fraction. 

Fig. 5: The position of (a) the maximum velocity and (b) the maximum 
temperature in the y — (f> plane with / = 100 and several values of A. 
The data are averaged over 200 cycles. Note that the two sets of the 
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curves are almost identical. 

Fig. 6: The maximum value of (a) the vertical velocity v y (y, (j) max (y)) and 
(b) the temperature T(y,(f) max (y)) at height y with / = 200. The data 
are averaged over 200 cycles. 
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